Временные ряды.
Неделя 3
Неделя 3
from google.colab import drive
drive.mount('/content/gdrive',force_remount=True)
cd /content/gdrive/'My Drive'/Colab_Notebooks/ML_c6_TAXI
PATH_TO_DATA = 'data'
import warnings
warnings.filterwarnings('ignore')
from sklearn.base import BaseEstimator, TransformerMixin
import os
import pandas as pd
import numpy as np
from scipy import stats
from tqdm import tqdm_notebook
from matplotlib import pyplot as plt
import seaborn as sns
sns.set_style("darkgrid")
import gc
import statsmodels.api as sm
from sklearn.linear_model import LinearRegression
from itertools import product
def prepare_dataset_geo_time_series(taxi_data,reg_data):
"""
Описание
Данная функция принимает на вход загруженные
таблицу с поездками и список зон с координатами
Возвращает временной ряд частот поездок по времени
размера (time_length х regions_length)
"""
df = taxi_data.copy()
regions = reg_data.copy()
#границы города, y - широта, x - долгота
(Ymin,Ymax) = (regions.south.min(),regions.north.max())
(Xmin,Xmax) = (regions.west.min(), regions.east.max())
#очистка
df.drop(df[df.pickup_latitude <= Ymin].index,inplace =True )
df.drop(df[df.pickup_latitude >= Ymax].index,inplace =True )
df.drop(df[df.pickup_longitude <= Xmin].index,inplace =True )
df.drop(df[df.pickup_longitude >= Xmax].index,inplace =True )
df.drop(df[df.passenger_count == 0].index,inplace =True)
df.drop(df[df.trip_distance == 0].index,inplace =True)
df.drop(df[df.tpep_pickup_datetime == df.tpep_dropoff_datetime].index,inplace =True)
#столбец времени по часам для построения ряда по часам
df['time'] = df.tpep_pickup_datetime.apply(lambda x: x[:-6])
#построение таблицы со строками - часами, столбцами - зонами,
#на пересечении - количество поездок в данный час из данной зоны
binx = sorted(list(regions.west.unique()) + [regions.east.max()])
biny = sorted(list(regions.south.unique()) + [regions.north.max()])
hour_counts = []
time = []
for (tt,sub_df) in tqdm_notebook(df.groupby('time')):
x = sub_df.pickup_longitude
y = sub_df.pickup_latitude
res = stats.binned_statistic_2d(x,y,x, statistic = 'count', bins = [binx,biny])
counts = res.statistic.reshape(-1)
time.append(pd.to_datetime(tt))
hour_counts.append(counts)
#датафрейм-результат, отсортированный по времени в формате datetime
time_series = pd.DataFrame(hour_counts,index = time,columns = list(regions.region),dtype=int)
time_series.sort_index(inplace=True)
return time_series
%%time
regions = pd.read_csv(os.path.join(PATH_TO_DATA,'regions.csv'),sep=';')
regions.head()
'''
%%time
months = ['2015-12','2016-01','2016-02','2016-03','2016-04','2016-05']
time_series_tables = []
for month in months:
df = pd.read_csv(os.path.join(PATH_TO_DATA,f'yellow_tripdata_{month}.csv'))
time_series_table = prepare_dataset_geo_time_series(df,regions)
gc.collect()
time_series_tables.append(time_series_table)
del df
time_series_6_month = pd.concat(time_series_tables)
#time_series_6_month.to_csv('time_series_6_month.csv')
'''
Загрузим данные за 6 месяцев
%%time
time_series_6_month = pd.read_csv(os.path.join(PATH_TO_DATA,'time_series_6_month.csv'),index_col=0)
%%time
time_series_6_month.columns = range(1,2501)
dt = [pd.to_datetime(tt) for tt in list(time_series_6_month.index)]
time_series_6_month.index = dt
time_series_6_month.head()
time_series_6_month.info()
Будем рассматривать зону, где расположено Empire State Building
ESB_lat = 40.748817
ESB_long = -73.985428
print(regions[regions.south<ESB_lat][regions.north>ESB_lat][regions.west < ESB_long][regions.east > ESB_long])
ESN_reg=1231
ESB_ts = time_series_6_month[ESN_reg]
with plt.xkcd():
plt.figure(figsize=(50,6))
plt.plot(ESB_ts, color = 'blue',linewidth = 1.5)
plt.grid(True)
plt.xlabel('time')
plt.ylabel('trip count')
Посмотрим на результаты stl декомпозиции для фрагмента ряда
def stl_plot(ts):
with plt.xkcd():
result = sm.tsa.seasonal_decompose(ts)
plt.figure(figsize=(50,6))
result.observed.plot(title = 'observed')
plt.grid(True)
plt.figure(figsize=(50,6))
result.trend.plot(title = 'trend')
plt.grid(True)
plt.figure(figsize=(50,6))
result.seasonal.plot(title = 'seasonal')
plt.grid(True)
plt.figure(figsize=(50,6))
result.resid.plot(title = 'residuals')
plt.grid(True)
stl_plot(ESB_ts[-500:])
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(ESB_ts)[1])
На графике тренда видна недельная сезонность, на графике сезонности видна суточная сезонность. Зададим регрессионные признаки, чтобы описать недельную сезонность
def reg_feats(start,stop,kw=2,ky=0,linear_fe=False):
#Вспомогательная функция для построения признаков
arg = np.array([x for x in range(start,stop+1)])
fe_list = []
columns = []
#Линейный признак для описания простейшего тренда
if linear_fe:
fe_list.append(arg)
columns.append('t')
#Признаки для описания недельной сезонности
for i in range(1,kw+1):
fe_list.append(np.sin(arg*2*np.pi*i/168))
columns.append(f'sin(2*pi*t*{i}/168)')
fe_list.append(np.cos(arg*2*np.pi*i/168))
columns.append(f'cos(2*pi*t*{i}/168)')
#Признаки для описания годовой сезонности
for j in range(1,ky+1):
fe_list.append(np.sin(arg*2*np.pi*j/8766))
columns.append(f'sin(2*pi*t*{j}/8766)')
fe_list.append(np.cos(arg*2*np.pi*i/8766))
columns.append(f'cos(2*pi*t*{j}/8766)')
df = pd.DataFrame(fe_list).T
df.columns = columns
return df
Возьмем количество регрессионных признаков 40 (kw=20) для описания недельной сезонности.
feats = reg_feats(1,ESB_ts.shape[0],kw=20,ky=5)
feats.head()
Построим линейную регрессию и будем анализировать остатки.
regressor = LinearRegression()
regressor.fit(feats,ESB_ts)
preds = regressor.predict(feats)
plt.figure(figsize=(50,16))
ax=plt.subplot(211)
ax.plot(ESB_ts.index,preds)
plt.title('regression');
resids = ESB_ts - preds
ax=plt.subplot(212)
ax.plot(resids)
plt.title('regression residuals');
ts_df = resids.to_frame(name="resids")
stl_plot(ts_df.resids[-500:])
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(resids)[1])
В остатках видны тренды и суточная сезонность, при этом недельная сезонность почти устранена
ts_df['diff1_resids'] = ts_df.resids - ts_df.resids.shift(24)
stl_plot(ts_df.diff1_resids[-500:])
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(ts_df.diff1_resids[25:] )[1])
ts_df['diff2_resids'] = ts_df.diff1_resids - ts_df.diff1_resids.shift(24)
stl_plot(ts_df.diff2_resids[-500:])
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(ts_df.diff2_resids[49:] )[1])
ts_df['diff3_resids'] = ts_df.diff2_resids - ts_df.diff2_resids.shift(24)
stl_plot(ts_df.diff3_resids[-500:])
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(ts_df.diff3_resids[73:] )[1])
ts_df['diff4_resids'] = ts_df.diff3_resids - ts_df.diff3_resids.shift(24)
stl_plot(ts_df.diff4_resids[-500:])
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(ts_df.diff4_resids[100:])[1])
ts_df['diff5_resids'] = ts_df.diff4_resids - ts_df.diff4_resids.shift(1)
stl_plot(ts_df.diff5_resids[-500:])
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(ts_df.diff5_resids[101:])[1])
Видно, что после 1 сезонного дифференцирования сезонная компонента в остатках регрессии имеет небольшую амплитуду (размах около 30). Остальные дифференцирования результат не улучшили, поскольку растет амплитуда сезонности. Можно, конечно, возразить, что уменьшается выраженность трендов, но данный факт можно проверить при дальнейшем подборе в рамках проекта. Пока остановимся на 1 сезонном дифференцировании.
plt.figure(figsize=(15,10))
ax = plt.subplot(211)
sm.graphics.tsa.plot_acf(ts_df.diff1_resids[25:].values.squeeze(), lags=168, ax=ax);
ax = plt.subplot(212)
sm.graphics.tsa.plot_pacf(ts_df.diff1_resids[25:].values.squeeze(), lags=168, ax=ax);
Возможные параметры Q=1,q=3,P=1,p=1(из-за вычислительных возможностей будем брать лаги с корреляцией около 0.4 для обычных компонент и 0.2 для сезонных)
ps = range(0, 2)
d=0
qs = range(0, 4)
Ps = range(0, 2)
D=1
Qs = range(0, 2)
parameters = product(ps, qs, Ps, Qs)
parameters_list = list(parameters)
len(parameters_list)
Для быстрого подбора модели напишем класс для моделей с интерфейсом fit predict
class TimeSeriesModel(BaseEstimator, TransformerMixin):
'''
Модель для предсказания временного ряда.
Использует регрессию на тригонометрические признаки и модель SARIMA
'''
def __init__(self,kw=2,ky=0,linear_fe=False,params=(1,1,1),seasonal_params=(1,1,1,24)):
self.kw = 2
self.ky = ky
self.linear_fe = linear_fe
self.params=params
self.seasonal_params=seasonal_params
def reg_feats(self,start,stop,kw=2,ky=0,linear_fe=False):
#Вспомогательная функция для построения признаков
arg = np.array([x for x in range(start,stop+1)])
fe_list = []
columns = []
#Линейный признак для описания простейшего тренда
if linear_fe:
fe_list.append(arg)
columns.append('t')
#Признаки для описания недельной сезонности
for i in range(1,kw+1):
fe_list.append(np.sin(arg*2*np.pi*i/168))
columns.append(f'sin(2*pi*t*{i}/168)')
fe_list.append(np.cos(arg*2*np.pi*i/168))
columns.append(f'cos(2*pi*t*{i}/168)')
#Признаки для описания годовой сезонности
for j in range(1,ky+1):
fe_list.append(np.sin(arg*2*np.pi*j/8766))
columns.append(f'sin(2*pi*t*{j}/8766)')
fe_list.append(np.cos(arg*2*np.pi*i/8766))
columns.append(f'cos(2*pi*t*{j}/8766)')
df = pd.DataFrame(fe_list).T
df.columns = columns
return df
def fit(self,ts):
self.ts = ts
#Признаки для регрессии
feats = self.reg_feats(start=1,stop=self.ts.shape[0],kw=self.kw, ky=self.ky,linear_fe=self.linear_fe)
self.regressor = LinearRegression()
self.regressor.fit(feats,self.ts)
preds = self.regressor.predict(feats)
#Получаем остатки регрессии и обучаем на них SARIMA
resids = self.ts - preds
self.model=sm.tsa.statespace.SARIMAX(resids.values, order=self.params,
seasonal_order=self.seasonal_params).fit(disp=-1)
return self
def predict(self,start,stop,verbose=False):
#Получаем признаки регрессии для предсказания
feats = self.reg_feats(start,stop,self.kw,self.ky,self.linear_fe)
#Предсказание регрессии
pred_regression = self.regressor.predict(feats)
#Предсказание SARIMA
pred_resids = self.model.predict(start=start, end=stop)
#Итоговое предсказание как суперпозиция
result = pred_resids + pred_regression
if verbose:
plt.plot(pred_regression,color='r') #график регрессии для отладки
plt.plot(pred_resids,color='g')#график предсказанных остатков для отладки
return result
Опробуем модель на простых параметрах
model = TimeSeriesModel(kw=2,ky=0,linear_fe=False,params=(0,0,1),seasonal_params=(1,1,1,24))
%%time
model.fit(ESB_ts)
prediction = model.predict(ESB_ts.shape[0]-500,ESB_ts.shape[0]-1)
plt.figure(figsize=(50,20))
plt.plot(ESB_ts[-500:].index,prediction,color='r')
plt.plot(ESB_ts[-500:].index,ESB_ts[-500:].values);
Можно увидеть, что даже простая модель описывает сложную сезонность в обучающих данных, при этом на тесте, естесственно такого не будет. Сравним результаты разработанной модели и предложенной SARIMAX с регрессионными признаками.
%%time
model2=sm.tsa.statespace.SARIMAX(ESB_ts.values, exog=model.reg_feats(1,ESB_ts.shape[0],kw=2,ky=0,linear_fe=False), order=(0, 0, 1),
seasonal_order=(1, 1, 1, 24)).fit(disp=-1)
prediction2 = model2.fittedvalues[-500:]
plt.figure(figsize=(50,20))
plt.plot(ESB_ts[-500:].index,prediction,color='r')
plt.plot(ESB_ts[-500:].index,prediction2,color='g');
Можно заместить, что даже при такой простой постановке написанный класс эффективней предложенной модели в 2 раза, при этом, как видно из последнего графика, предсказанные величины практически идентичны. Поэтому будем далее использовать разработанную модель. Также для простоты при подборе будет рассмотривать не весь ряд, а только 720 последних часов.
%%time
results = []
best_aic = float("inf")
warnings.filterwarnings('ignore')
for param in tqdm_notebook(parameters_list):
#try except нужен, потому что на некоторых наборах параметров модель не обучается
try:
model = TimeSeriesModel(kw=20,ky=0,linear_fe=False,params=(param[0], d, param[1]),seasonal_params=(param[2], D, param[3], 24))
model.fit(ESB_ts[-720:])
#выводим параметры, на которых модель не обучается и переходим к следующему набору
except ValueError:
print('wrong parameters:', param)
continue
aic = model.model.aic
#сохраняем лучшую модель, aic, параметры
if aic < best_aic:
best_model = model
best_aic = aic
best_param = param
results.append([param, model.model.aic])
warnings.filterwarnings('default')
result_table = pd.DataFrame(results)
result_table.columns = ['parameters', 'aic']
print(result_table.sort_values(by = 'aic', ascending=True).head())
Обучим модель на данных за 5 месяцев, при этом последние 720 часов используем для валидации
model = TimeSeriesModel(kw=20,ky=0,linear_fe=False,params=(1,0,3),seasonal_params=(1,1,1,24))
%%time
model.fit(ESB_ts[:-720])
print(model.model.summary())
Посмотрим на точность на обучающих данных (последние 720 часов обучения)
prediction0=model.predict(ESB_ts.shape[0]-1440,ESB_ts.shape[0]-721)
plt.figure(figsize=(50,20))
plt.plot(ESB_ts[-1440:-720].index,prediction0,color='r')
plt.plot(ESB_ts[-1440:-720].index,ESB_ts[-1440:-720].values);
Видно, что в целом модель описывает данные, при этом есть сложности с описанием отдельных пиков и впадин, которые, скорее всего, можно объяснить аномалиями.
Посмотрим на остатки SARIMA
plt.figure(figsize=(15,8))
ax = plt.subplot(211)
ax.plot(model.model.resid[-720:])
plt.ylabel(u'Residuals')
ax = plt.subplot(212)
sm.graphics.tsa.plot_acf(model.model.resid[25:].squeeze(), lags=168, ax=ax)
print("Критерий Стьюдента: p=%f" % stats.ttest_1samp(model.model.resid[25:], 0)[1])
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(model.model.resid[25:])[1])
Остатки напоминают визуально белый шум. Уровень значимости для теста Стьюдента имеет пограничное значение, можно считать, что в первом приближении остатки несмещены. Критерий Льюнга-Бокса и коррелограмма говорят о слабой автокорреливанности ряда. Критерий Дики-Фуллера отвергает гипотезу о нестационарности, но, в связи с большим количеством данных, стоит посмотреть на stl декомпозицию ряда остатков. При этом нужно учитывать, что при больших данных аппарат проверки гипотез склонен отклонять нулевые гипотезы в пользу альтернатив.
index=[pd.to_datetime(tt) for tt in ESB_ts.index[-500:]]
ts_resids= pd.Series(model.model.resid[-500:],index=index)
stl_plot(ts_resids)
Можно заметить, что ряд остатков скорее нестационарен в связи с тем, что небольшие тренды и суточная сезонность осталась, хотя ряд уже больше напоминает что-то хаотическое визуально. С этим, скорее всего, связана слабая автокоррелированность ряда. Посмотрим, как модель ведет себя на валидации.
prediction = model.predict(ESB_ts.shape[0]-719,ESB_ts.shape[0])
plt.figure(figsize=(50,20))
plt.plot(ESB_ts[-720:].index,prediction,color='r')
plt.plot(ESB_ts[-720:].index,ESB_ts[-720:].values);
Точность на данных за май несколько хуже, чем на обучении. Нужно также учитывать, что обучение происходило на данных за 5 месяцев, при этом не учитываются особенности месяцев в году, то есть годовая сезонность. В качестве улучшения, можно попробовать итерационно описывать остатки данного шага снова моделью ARIMA и просто добавлять результат при предсказании.